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We study theoretically Doppler laser-cooling of a cluster of 2-level atoms confined in a linear 
ion trap. Using several consecutive steps of averaging we derive, from the full quantum mechanical 
master equation, an equation for the total mechanical energy of the one dimensional crystal, defined 
on a coarse-grained energy scale whose grid size is smaller than the linewidth of the electronic 
transition. This equation describes the cooling dynamics for an arbitrary number of ions and in 
the quantum regime. We discuss the validity of the ergodic assumption (i.e. that the phase space 
distribution is only a function of energy). From our equation we derive the semiclassical limit (i.e. 
when the mechanical motion can be treated classically) and the Lamb-Dicke limit (i.e. when the 
size of the mechanical wave function is much smaller than the laser wavelength) . We find a Fokker- 
Planck equation for the total mechanical energy of the system, whose solution is in agreement with 
previous analytical calculations which were based on different assumptions and valid only in their 
specific regimes. Finally, in the classical limit we derive an analytic expression for the average 
coupling, by light scattering, between motional states at different energies. 
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I. INTRODUCTION 



The development of laser cooling and trapping techniques in the last decades has allowed for many spectacular 
experimental achievements jjj. Hand in hand with the experimental work, a theory of laser cooling has been widely 
developed, providing a precise description of many experimental situations (^||. Yet, the great majority of the 
theoretical studies have dealt with the interaction of laser light with single particles (being the same as many non- 
interacting particles). The treatment of a many-body system coupled to light is considerably more complex, due to 
the large number of degrees of freedom which results from the interaction among the constituents of the system. 

In this work we investigate laser cooling of a many-body system, taking as a representative - and experimentally 
relevant - example a one-dimensional Coulomb cluster, i.e. a crystallized structure of ions which are trapped in a 
Paul or Penning trap We develop a model for the dynamical behaviour of the crystal's mechanical energy. In 

this system the physical processes in play are the trapping potential which confines the ions, their mutual Coulomb 
interaction, and their interaction with laser light. For sufficiently small kinetic energy the motion of the ions is 
properly described by the collective excitations of the cluster, the so-called harmonic modes of the crystal, due to the 
interplay of trapping potential and Coulomb repulsion. The problem is then characterized by three main frequency 
scales: (i) the oscillation frequencies of these modes, which are determined by the trap frequency v and the number of 
particle N (for ions of equal mass and charge); (ii) the recoil frequency u>r characterizing the exchange of mechanical 
energy between radiation and atoms, lor = 2it 2 Ti / m\ 2 , with m mass of the ions and A wavelength of the light coupling 
quasi-resonantly to the electronic transition, and (iii) the linewidth of the electronic transition 7, which characterizes 
the rate at which photons are scattered by the atoms. Doppler cooling of particles in a harmonic trap focuses on the 
regime 7 > v. 

The theory is well developed for single trapped particles in specific regimes: in the semiclassical limit 7 3> v,ur, 
when the mechanical motion can be treated classically, and in the Lamb-Dicke limit v 3> lor, when the size of the 
motional wave function is much smaller than the laser wavelength || . In the semiclassical case, the standard procedure 
consists in writing the equations for the system dynamics in the Wigner representation, adiabatically eliminating the 
excited electronic state and then expanding in h up to the second order. In this way a Fokker-Planck equation in the 
position and the velocity of the ion is derived |j. In the Lamb-Dicke limit, corresponding to an expansion in ujr/u 
up to the first order, the cooling dynamics are described by a set of rate equations projected on the electronic ground 
state and on the eigenstates of the motion, which can be analytically handled ■ Both treatments yield in the limit 
7>i/ the same final distribution of cooled atoms in energy space, and the same dependence of the final temperature 
on the cooling parameters. Our approach shows how they are special cases of the same equation that we derive for 
the mechanical energy of the crystal. 



1 



The presence of many mutually interacting particles complicates the treatment considerably, because the mechanical 
Hamiltonian, having an increased number of degrees of freedom, often does not allow for a simple and transparent 
solution. An immediately visible effect arises in the spectrum of the mechanical energy, where the increased number 
of degrees of freedom is in general connected with the appearance of quasi-degeneracies and with a dense distribution 
of the energy levels pToj ] . The situation is facilitated again if one restricts the treatment to the Lamb-Dicke regime 
PUP^f but in the more general case we want to discuss, different steps of simplification can be made as will be shown. 

The treatment in this paper focuses on Doppler cooling of an iV-body one-dimensional Coulomb crystal, in particular 
on the type of linear ion chains obtained in linear Paul traps fi|l3|]. Here, the many-body mechanical Hamiltonian can 
be approximated by a set of harmonic oscillator modes, all having different frequencies u%, ...,i>n, where each mode 
corresponds to a collective excitation of the crystal. Since for Doppler cooling the rate of photon scattering, and thus 
the cooling dynamics is determined by 7, we define a grid in energy space AE <C ^7 and study the cooling dynamics 
among the energy shells defined by the grid. From this starting point, with a procedure of successive averaging we 
derive a rate equation for the population at each motional energy on the grid. From this result, the semiclassical 
limit is recovered by expanding in lur/^. In particular, we get a Fokker-Planck equation for the motional energy 
describing cooling of N ions, whose solution agrees with the ones of standard semiclassical treatments P, |l2] , |l4| , |l5f . 
In the Lamb-Dicke regime, starting from the full master equation, we discuss the conditions under which a compact 
equation for cooling can be found. Finally, we derive an explicit form of the rate equation in the limit 7 3> v by 
evaluating an analytic expression for the average coupling by light scattering between motional states at different 
energies. 

Although we restrict the investigation to a one-dimensional Coulomb crystal, the theory can be directly extended 
to three dimensions, and the results are applicable to clusters in both Paul and Penning traps. More in general, the 
idea behind our treatment is that - in incoherent processes - if the rate determining the dynamics of interest of the 
system can be singled out, the contribution of processes occurring on faster time scales is often well represented by 
the average value on the slower time scale, characterizing the incoherent process. 

This work is organized as follows. In Section II the model from which our derivation starts is introduced and discussed. 
In Section III we discuss and apply a series of approximation, from which we obtain an equation for the cluster's 
mechanical energy describing cooling of the crystal. In Section IV we study the equation in the Lamb-Dicke and 
Semiclassical limits, and we compare our results with those obtained in the same limits with other treatments. In 
Section V we derive an explicit functional form of the rate equation in the limit 7 S> v. Finally, we draw our 
conclusions. 

II. MODEL 

In this section we introduce the model which we study throughout the paper. We first discuss the mechanical properties 
of a one-dimensional Coulomb crystal, then the interaction of laser light with the internal electronic transition of each 
ion, and finally the mechanical action of the light on the collective mechanical degrees of freedom of the crystal. 

1. Mechanical properties 

The mechanical potential on an ion cloud in a Paul trap is the sum of the potential exerted by the trap on each ion 
and of the Coulomb repulsion among the ions. Sufficiently far away from the trap electrodes, the potential of a Paul 
trap can be considered harmonic and the total potential for N ions of charge e has the form: 
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where fj = (xj,yj, Zj) is the position of ion j and u J 0x , u 3 0y , u J 0z depend on the trap parameters and on the mass of 

the jth ion. At sufficiently low temperatures the ions crystallize around the classical equilibrium positions which 
are solutions of the set of equations dV/dfi]^) — ||. Sufficiently below the crystallization temperature the motion 

of the ions around these equilibrium points is harmonic to good approximation. In this limit, the total mechanical 
potential can be described by its Taylor series truncated at the second order in the expansion around {r^}. In a Paul 
trap with cylindrical geometry and very steep potential in the radial (x and if) direction the ions crystallize along 
the trap (z-) axis. In this regime the amplitude of the radial oscillations is much smaller than the axial ones. Here, 
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we assume the radial degrees of freedom to be frozen out such that the motion is one-dimensional along the z-axis, 
and we take the axial potential to be electrostatic, as it is the case in a linear Paul trap |ll|. Then the truncated 
mechanical potential has the form: 

f N 

V(q u ...,q N ) = 2 zJ V ^QjQk, (2) 

j,k=l 

where Vjk is a real, symmetric and non-negative matrix (the explicit form of its elements can be found, for example, 
in fill), while qj = Zj — z^ are the displacements from the equilibrium positions Given m mass of each ion, the 
secular equation for the harmonic motion of the crystal has the form: 

Y^Vijbf = mulbf with a=l,...,N (3) 

j 

where v a are the eigenvalues and 6" the associated eigenvectors, which are complete and orthonormal. 
In the basis q' a = bf qi the motion is described by the Lagrangian for N independent harmonic oscillators of 
frequency v a . Given the canonical momentum p' a = q' a conjugate to q' a , the motion is quantized by associating a 
quantum mechanical oscillator with each mode. Denoting by a a and aj, the annihilation and creation operators for 
the mode a, respectively, the coordinates qi are now written as 



«-E(Or^(-+4). («> 

a 

and the Hamiltonian for the mechanical motion has the form: 

H mec = ^ hv a {a) a a a + - J (5) 

The energy eigenstates of each mode a are the number states \n a ) with eigenvalues (energies) £„ Q = Ti(n a + ^)v a . In 
the following we use states |n) = \n 1: n 2: ...njv) to describe the eigenstates of the motion of the crystal corresponding 
to the eigenvalues £ n = J2 a ^n a ■ 

In summary, N trapped ions crystallized in one dimension can be described by N harmonic oscillators with frequen- 
cies fx, vn- These frequencies are solutions of (||) and, what is important, they are incommensurate. Therefore the 
spectrum of mechanical energies of the system does not exhibit the discreteness and equispacing property of the single 
harmonic oscillator spectrum, but shows a dense distribution of levels, and at sufficiently high energies and large N 
it assumes a quasi-continuum character. In this limit we can define the density of states of the system, g(E), which 
is a smooth function of the energy and is defined on a grid of energies AE, such that the number of states D(E) 
contained in the interval T{E) = [E — AE/2,E + AE/2] is D(E) > 1, and D(E) ~ g(E)AE. For N modes (ions) 
g{E) can be evaluated by solving the following integral: 
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E N-1 



E pEi pEw-i 



(nvi)(hv2)...[hv N ) AE J J 



n 



(6) 



i.e. by taking the derivative of the ratio between the volume in phase space of energy < E and the volume occupied 
by a single state, {hvi){hv2)---{JiVN)- In Fig. 1 we plot the number of states as a function of the total mechanical 
energy for a chain of three ions, taking a grid AE — hvi/5, and compare the exact value with the smoothed function 
g(E)AE, with g(E) given by (|). 
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FIG. 1. Number of states D(E) as a function of the total energy E/Tiv evaluated for three ions on a grid AE — hv/5 (dots 
joined by the grey line); smoothed function g(E)AE of Eq. ^ with N = 3 (black line). The frequencies of the modes are 
vi =v,V2 = 1.7321;/, v z = 2.4083;/. 



2. Interaction with light 

We consider that laser light with frequency u>l and wave vector k = (k x , k y ,k z ) drives some ions of the crystal, coupling 
to their internal two- level transition with electronic ground state |g), excited state |e) and resonance frequency uio. In 
the Rotating Wave Approximation and in the frame rotating with the laser frequency, the total Hamiltonian has the 
form: 



H = H a 



H 



AL- 



Here H a t is the Hamiltionian for the internal degrees of freedom, defined as 



N 

H at = -hSj2\e)j(e\ 



(7) 



(8) 



where S = wj, —ujq is the detuning of the laser from the atomic resonance and j labels the ions. H mcc is the mechanical 



Hamiltonian defined in (|5|), and Hal describes the interaction between laser and atoms, 



2 1 J 



(9) 



where fl(zj) is the Rabi frequency at the position Zj, cr~j~ , o~ are the raising and lowering dipole operators, respectively, 
for the internal state of the j th ion, while {j} is the set of driven ions. Assuming that the light intensity does not 
vary rapidly in the vicinity of z^°\ we can approximate £l(jZj) ~ ^(z^) — The operator e -lfe * Zj in (||) represents 
the mechanical effect of the light interaction, i.e. a shift in momentum space by one photon recoil which goes along 
with the excitation. Using (|J), its explicit form is: 



exp(— ik z Zj) — exp —ik z z 



(0) 



N 



(10) 



where rjf is the Lamb-Dicke parameter for mode a and ion j, defined as p 



2mv a 



(11) 



For some ions we can have flj — 0, which means that not all ions of the chain are driven. Such condition can be 
achieved either when the ions are sufficiently spaced to allow for their individual addressing Jl8[ ] , or by introducing 
a different type of ion into the crystal, such as an isotope of the trapped element |l9|, or a different species |L7 20 
whose transition frequency is not resonant with the laser. In this case, some differences in the mechanical behaviour 
arise, which, however, are not relevant for the results that we derive below. 
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3. Master Equation 



We describe the dynamics of the driven crystal through the master equation for the density matrix p of the iV-ion 
system: 

±p = -L[H,p]+L(p), (12) 

where H has been defined in (Q) , and L is the Liouvillian describing the incoherent evolution of the system due to its 
coupling to the modes of the electromagnetic field: 

AT 

Hp) = i », 4*7 p 9*1 ( 13 ) 

Here 7 is the decay rate of the internal excited states \e)j, and pj is the density matrix describing the feeding for the 
ions j, shifted in momentum space by the recoil of the spontaneously emitted photon, 

p 3 = j d(cos6)Af(cos8)e lkcosez >pe- lkcos9z >, (14) 
with AT (cos 8) being the dipole pattern of the spontaneous decay, and k = \k\. 



4- Low saturation limit 



In the limit of low saturation, Q <C 7, the excited states \e)j can be eliminated from ([12]) in second order perturbation 
theory pl| ]. Thereby we obtain a closed equation for the internal ground state |g) = Oj=i \&)j> wn ^ cn we project on 
the basis of mechanical states |n). The equations describing the evolution of this system are |10|: 



d i 

— (n|p|m) = -^(^n - £ m )(n\p\m) 



(15) 



^ r(n|e^^|k)(k|e- l ^^|l) . . . , (l|e ifc *^|k>(k|e- ifc *^|m) 
— ~ ; ^-rr(l\p\m) 



M k,l 
•>2 



(£ k - £x)/h - 6 - if/2 



(f k - £i)/h - S + if/2 



n|p|l> 



EjE f d(cos^(cos0)(n 



"-'|k)(k|e- ifca ^|r)(s|e ifca ^ |j)(j|e- !;fccos ^|m)(r|p|s) 



Ak cos 9zj 



[(£j - £ s - £ k + £ r )/fi + i 7 ][(£j - £„)/fi - S - i^y/2] [(£ k -£ r -£ s + £ s )/h - if][(£ k - £ T )/h -5 + if/2]_ 

Here, the coefficients describing the coupling among the populations, (n|p|n), and the coherences, (n|p|m), are 
proportional to the Franck-Condon coefficients, which have the form: 

N 



(l|e-^|n) = n(U e ~ !;<(aQ+4) K> 

a=l 
N 

= n ° xp 



nr 

2 



(r a + \l a - n a \)\ 



(16) 



where r a = min(Z a , n a ), and lJ-'° (x) is a generalized Laguerre polynomial. They represent the probability 
amplitude of a transition from the initial motional state |n) to the final motional state |1) by absorption or emission 
of a photon. 
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III. AN ERGODIC EQUATION FOR LASER COOLING OF THE CRYSTAL 



In the following, starting from the full quantum mechanical Eqs. ( |15| ) and using consecutive steps of averaging, we 
derive an equation for Doppler cooling of an ion crystal. This equation will be ergodic, in the sense that it describes 
the dynamics of the system in terms of the population P(E) at the crystal's mechanical energy E. 

In the Doppler cooling limit (which implies 7 > v a for all a = 1, N, and <5 = 0(7)) we define a grid of energies 
AE such that in an interval T{E) — [E — AE/2,E + AE/2] the number of states D{E) ^> 1, and the Lorentzian 
describing the resonant response of the atom varies infinitesimally on each interval. The resulting shells of energies 
can be visualized as shown in Fig. 2 for the case of a two-ion crystal, where each axis represents the energy of one 
mode, and each point is a state |n) = r^a) . Here, the states with total energy £ fall into the shell of width AE 
along the line E = Kv\{rii + |) + ^2(^2 + j)- 



E' 



E 




FIG. 2. Coarse grained energy space for the case of two modes of frequency v\, U2- The points are the states with energy 
£n = £nx + £n 2 - The broad lines represent two energy shells J-(E) and J-(E'). 



We can now decompose the sums over the motional states appearing in (g_5|) as follows: 

E=EE' where E = E 

k E k k k k : £ k G^(-E k ) 



(17) 



where now E^ is one energy value on the grid, and J-"(i?k) = [Ek — AE/2, E^ + AE/2] is the energy interval centered 
in Ek with width AE. Using these definitions we rewrite the equations for the populations in ( |l5| ) and sum over all 
states that lie within the same energy shell JF(i?k): 



at 



4 5[(^k-Sn)A-<5] 2 +7 2 /4 n , k 

-E?E 

{3} E k ,E, 
>2 



1 



(Ek - E\)/h - 5 - i 1 /2 



£'|(n|e^|k)| 2 (nHn 

n,k 

£'£'<n|e ifc ^|k><k|e 



— ik z 



|l)(l|p|n) 



n,k 1,1/n 



+7^7 / d(cos0)AA(cos0) Y, 



{3} 

EE' E 

-Bj,B s n,k,s 



[(£ k - £ r )A - c + 7 2 /4 



ik cos Qzj 



k)n<k|e-^|r)|>Mr) 



ne 



ifc cos 6z 



J |k) (klc 



r)(s|e^^|j)(j|c- 



--iA; cos Oz 



' n) 



[{E s -E s -E k + Er)/n + i~f] [{Ej -E s )/n-S- i 1 /2\ 



r|p|s 



j^k or r^s 



(18) 

(19) 
(20) 
(21) 



where we have separated the terms involving the populations, ( [iq ) and (p0|), from the ones which involve the co- 
herences, ( p^| ) and (^l|). We now compare the term (18), the loss rate of the population, with the coupling to the 
coherences in (|l9|): Line ( |T8| ) is proportional to the sum of moduli square of Franck-Condon-coefficients, thus it is 
the sum of positive terms. The sum in line ( |T9| ) adds up coefficients with alternating signs, as can be verified from 
Eq. (|l6|). For D(E) > 1 we can therefore safely assume that 



G 



^^'^'|(n|e^^|k)| 2 (n|p|n)»|^ £ '(n|e^ |k) (k|e"^ |1) (l|p|n) | (22) 

{j} n k {j} n,k,l#n 

This is equivalent to a Random Phase Approximation. On the basis of this consideration the coupling between 
populations and coherences in (|l^) may be neglected in comparison to the loss rate described by the term in (|Tg|). 
Analogously we can neglect the terms in ( ^lf ) in comparison with (po|), and we obtain the following set of rate equations: 

Ay'(nHnH-7?yy " 3 T '(n|p|n)V '| (n|e 4 *«* |k) | 2 (23) 

d * n 4 [(^k-Sn)A-5] 2 +7 2 /4^ k 



4 m- 7 - 1 £^ [(^k - # r )A - ^] +7-/4 " k 



where all Rabi frequencies have been set equal, Qj = Q. This simplifies the discussion in the next subsection, but 
the more general case of different Qj can be treated with the same methods under conditions which will be discussed 
below. 



A. Ergodic hypothesis 

Equation (^) describes how laser light couples to the internal electronic transition and exchanges mechanical energy 
with the crystal. The probability amplitude for the radiative event to occur is weighted by a Lorentzian distribution, 
which below saturation has width 7, and which is a function of the difference between the mechanical energies of the 
crystal before and after the scattering. The mechanical effect of light on the motion of the crystal is described by the 
operator (|Io|), and the probability for the crystal initially in state |n) to be scattered into state |k) by the absorption 
or emission of a laser photon is given by the modulus square of the Franck-Condon coefficient, Eq. (|l6|). When we 
regard the modulus square of Eq. ( |l6|) as a distribution over the states |k) after an absorption or emission event, given 
the initial state |n), then the average motional energy transferred to the ion cluster and its variance are calculated as 
(see Appendix A): 

(£k - £„)k = E( £ k - £n) \(n\e* kcos9z i\k)\ 2 = hu, R cos 2 6, (24) 

k 

Of 

({£ k - £ n) 2 )k = hu R cos 2 6^- + (Hlur cos 2 6) 2 , (25) 

where lor = hk 2 /2m is the recoil energy. Thus, the first and second moments of this distribution depend on the 
(single-ion) recoil energy lor and on the energy per ion of the initial state, but not on the details of the quantum state 
In). 

Guided by this result, we will now make the approximation that in the limit of Dopplcr cooling, 7 > Vj, and for 
sufficiently large density of states, D(E) ^> 1, the oscillations of the Franck-Condon coefficients with the vibrational 
numbers play a negligible role, whereas the average properties determine the cooling dynamics. We formulate this 
assumption by first introducing an average Franck-Condon coupling between shells of energy E n and E^, 

1 N 



and then writing 



1 N 

-EE'Kkl^l^l^^^Q^^n^k), (27) 
j=l k 

thereby neglecting the dependence of the LHS on the details of the state |n). We will discuss this assumption and 
possible alternatives in the next section. It will be one of the important results later to determine an explicit expression 
for Q (k) {E n , Ek) in Eq. @. Furthermore, in @ all ions of the crystal are assumed to be driven. In the next section 
we will discuss the case in which only a set is driven. 

Applying (p7|) to Eq. (E3), in the limit in which all ions are driven, we obtain 
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1 V' (n|p|n) = -N^ V VVlpln) 

+ 7V 7 / d(cos0)AA(cos0) ^ — 2 2^ r H r > 

4 E t^ r [{Ek- E r )/h- 5} + 7 2 /4 



(28) 



In order to obtain the second line of ( pq ) from the second line of ( p3| ) we have made a further approximation: First 
we write 

^|(k|e-^^|r)| 2 |(n|e lfccos ^|k)| 2 = ^ | (k| e -^ 2 ^ |r) | 2 ^ | (n|e lfc cosSzi |k)| 2 6 jU (29) 
j j i 

and then we set Sji — > 1/iV. This corresponds to the assumption that, in the regime we are considering, the mechanical 
effect of the process of absorption + emission of a photon by one ion is equivalent to absorption by one ion and emission 
by another ion, weighted by the probability for the two ions to be the same. 
By defining the population densities P(E, t) of the energy shells through 

AEP(E ni t)=J2'Hp\n) , (30) 

11 

such that J °° dEP(E, t) = 1, and by using D(E) = g(E)AE, we finally arrive at a rate equation as a function of the 
motional energy: 



—P(E,t)=-j—N dE 1 yy v ' iy 

df 4 Jo [( Sl _ S )/^_j] 2 +7 2/4 



-P[E,t)=- 1 — N\ dE l r/ J y ^ L» lL— P(E,t) (31) 



+ ,^N f 1 d(cos W cos,) r^.^Q^-)^,^^^^)^,^) 
4 7_! Vo Jo {(E 1 ~E 2 )/h-5} 2 + 7 2 /4 

where we have replaced the sums over the energy by integrals, valid when the average energy (E(t)) 3> AE. The two 
parts of the rate equation (|3l| ) show how, after the various steps of averaging, the total population of a shell with 
energy E changes in time: Population is lost by excitation to shells with energy E\ and subsequent emission into any 
shell, and population originally at Ei is excited to a shell at E\ and then scattered into a shell at energy E. 

It should be noted that the restriction to a one-dimensional crystal enters into the number of modes and into the 
geometry of the laser beams, as well as into the pattern of the emitted radiation. However, the application of the 
"Random Phase Approximation" (|22|) and of assumption (|27]) , which are at the basis of our derivation, are in no way 
restricted by the dimensionality of the problem. Actually, an extension to three dimensions would endorse the two 
approximations, since the number of modes increases, and with it the density of states in the spectrum of the motional 
energies. Therefore, an equation for the total motional energy of the crystal of the form of |Il]) can be derived in 
three dimensions using the same considerations done here for the one-dimensional case. 

B. Discussion 

1. Ergodic assumption 

Since the distribution of the total mechanical energy over the system is sufficient for describing the cooling process, 
assumption (p7|), leading to Eq. ([H]), simplifies the description of the laser cooling dynamics of the crystal significantly 
by reducing dramatically the dimensionality of the problem. Equation ( |3l| ) could have been also obtained by assuming 
that the population of all states in the same energy shell is equal, 



(n\p\n) = p(E n ), (32) 



thus defining p(E) as the population of a state at energy E, such that (c.f. Eq. (30)) 



P(E)=p(E)g(E). (33) 



8 



Equation ( |32| ) implies ( p7| ) and is an even stronger assumption. It corresponds to assuming that the system behaves 
ergodically, i.e. that the populations of the states in an energy shell equalize faster than the average quantities of the 
system evolve in time. It leads to the same rate equation ([3l]), but it is a more natural assumption in the last stages 
of the evolution, when the system tends asymptotically to the steady state. 

In general, ( |3l| ) describes the dynamics of cooling to the extent that assumption ( p7| ) is valid, i.e. when the average 
coupling between an eigenstate of the mechanical energy and the states of an energy shell is a smooth function of their 
respective energy. This is true when the average coupling of each state of one energy shell to the states of another 
energy shell is of the same order of magnitude. Outside the Lamb-Dicke regime, given the oscillatory behaviour of the 
Franck-Condon coefficients, this holds if each state couples appreciably to more than one state of the other energy 
shell, i.e. for D(E) ^S> 1 and 7 sufficiently large. In the Lamb-Dicke regime analogous considerations can be applied, 
and they will be discussed in detail in the next section. 

The ergodic regime can also be justified by the existence of a physical process whose main effect is to thermalize 
the states within one energy shell and which acts on a shorter time scale than the energy-changing processes. This 
assumption is at the basis of treatments in the kinetic theory of quantum gases [^2| , where the interatomic collisions 
lead to thermalization, and the gas can be considered in a thermal quasi-equilibrium distribution on the time scale 
in which it is cooled. The assumption of rapid thermalization is also central to an earlier study of laser cooling of 
Coulomb clusters There the effect of mode- mode coupling (anharmonicity) has been proposed as a possible 

agent, which does not explicitly appear in the equations but justifies condition (\Z7\). Direct evidence of this effect has 
been found in numerical studies in |17| for the case of exact degeneracy between the modes frequencies. 



2. On the laser intensity distribution over the crystal 

In deriving ( |3l| ) we have assumed that all ions are uniformly driven. We discuss now the case in which only a subset 
{j} of ions in the crystal are illuminated. For simplicity, let us assume in Eq. (|3|) that fi; = fi for I £ {j}, and fli = 
otherwise. Due to the geometry of the crystal, the coupling of each ion to a mode is a function of the mode and of 
the ion's position in the crystal [see Eq. (JO)]. Some positions in the chain are more strongly coupled to a certain 
set of modes, and some positions are even decoupled from certain modes [ p"7| . Then, the ergodic assumption is valid 
provided a suitable set of ions is driven, such that the coupling of each mode of collective motion to the radiation is 
of the same order of magnitude. In that case, the ergodic assumption, Eq. (|27|), can be expressed as: 

i^^'|(k|e-^|n}| 2 ^(B k )Q W (£n,B k ), (34) 

« k 

where M is the number of driven ions (M < N) . Taking into account this more general case, the rate equation ( |3~l| ) 
has the form: 

giEjQM&E!) 



dt y '4 J [{E 1 -E)/h-5f + 7 2 /4 V ' 



7 ?M t d(cos W cos,) /°V /^ £ /^ mfl) (^%^f fe) (^^) p(£2;t) , (35) 



1 Jo Jo ' [{E l -E 2 )/h-8Y +7 2 /4 



Finally, note that the number of driven ions, M, appears in (J3lJ) and (|35|) as an overall scaling factor, making the 
optical pumping rate of the crystal M times that of a single ion. In fact, it has been shown earlier that below 
saturation the ions behave as independent scatterers, and their contributions to cooling simply add up [ ^0|Jl7| . 



IV. THE SEMICLASSICAL LIMIT AND THE LAMB-DICKE LIMIT 

In this section, starting from Eq. (^lj) or ( |35| ) we derive the limit where the motion can be treated classically (semi- 
classical limit), and obtain a Fokker-Planck-cquation for the energy which is analytically solvable. We compare our 
result with the well-known treatments of || and and find full agreement with their theoretical predictions. In 
the second part of this section, starting from (15) we consider the Lamb-Dicke limit, and discuss the conditions under 
which an equation for the motional energy of the type of ( pf| ) can be derived for describing cooling of N ions. 
For clarity in the derivation, we rewr ite Eq. (|l]) as: 
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j poo 

g(E)-p(E,t) = -g(E)p(E,t) J dE 1 f^ ) (E 1 )L(E 1 - E) 

/l p CO pOO 

d(cos^(cos9) / d£?i / dE 2 f^ cos9 \E 1 )f^\E 2 )L(E 1 -E 2 ) P (E 2 ,t). (36) 
-1 Jo Jo 



Here, p(-E) is the population of a state of energy E defined in (|32|), and L(£7i — i?) is the Lorentzian distribution, 



L(Ei — E) = 

A[{E l -E)/h-5Y + 7 2 ' 

while f E k \E') is defined through: 

ff ( J B')Q (fc) (^,i?') = /i fe) (^)- 

From this definition, using ( p7|) and ( p4f]25| ) , it can be verified that fg{E') fulfills the relations: 

dE'fg cose) (E') = 1, 



} 2E 
N 



pOO 

/ dE'(E' - E)f { E kcos0 \E r ) = hcoRCOS 2 9, 
Jo 

/>oo 

/ dE'(E' - E) 2 f { E kcos8) (E' 
Jo 



) = (Tilur cos 9) + hujR cos 



(37) 

(38) 

(39) 
(40) 
(41) 



The analogous relations follow for f^ z \E') when 9 is replaced by 9$, where k z — fccos#o- 
Equations ( [io|) and (^l|) are equivalent to (|24|) and ( p5| ) in the approximation of a smooth energy scale. 



A. Derivation of a Fokker-Planck equation for the energy 

We first consider the limit in which the recoil frequency is much smaller than the linewidth, ujr <C 7. In this limit. 



while fjp^E') varies on the scale hujR (see Eqs. (40-4l|)), the variation of the population p(E) on the energy scale 
Tiujr is small, i.e. 



hw R \-^p(E,t)\ <&p(E,t). 



We now expand Eq. (36) around E up to first order in the parameter u>r/j and, using (]3l^[l"l"|), we obtain: 
d 



dt 



g{E)p{E 1 t) = -2^cos 2 o L'(O)g(E)p(E,t) 



(42) 



(43) 



E 

{a + cos 2 9 )L(0) - 2— cos 2 6 L'(0) 



g(E)p'(E,t) 



w fi -|(a + cos 2 9 )L(0)g(E)p"(E, t). 



Here. 



and 



1/(0) = 



±L{x) 
dx 



4A/r2 2 7<5 



d^os^cos 2 ^^^) . 



(44) 



(45) 



We rescale the time as 



t = wr(cos 2 9 + a)L(0) t 



(46) 
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and define 

c 2 C ° s2g ° L/(0) f471 
(cos 2 O + a) L(0)' 1 ' 

Note that that with definition (|44|), C < for red detunings <5 < (w^ < ujq). Eq. ( [l3| ) can now be written as: 

^g(E)p(E, t)=~ [Ng(E)p(E, r) + Eg(E)p'(E, r)} + 1 [Ng(E)p'(E, r) + Eg(E)p"(E, r)} . (48) 

Using the energy dependence of the smoothed density of states g(E) as given in and (|33|), we get a Fokker-Planck 
Equation of the form: 

£pP{E,T) = ~ [A(E)P(E,t)} + \^ [B(E)P(E,r)] , (49) 

where 

ME) = 1 + £e, 
B{E) = 1e. 

Eq. (^) can be easily solved. In the following, we calculate the steady state solution and the time evolution of the 
system, and compare the result with the existing treatments evaluated in the limit of one ion (N = 1). 



1. Steady State 

The steady state distribution Pq(E) satisfies the equation -^P (E) = 0. Thus, it is solution of the differential equation: 

A [B(E)P (E)] = 2A(E)P (E), (50) 

(where the integration constant has been set to for a convergent solution) and has the form 

P (E) = FE N ~ 1 exp (CE) , (51) 

with normalization constant F. For red detunings (C < 0) and if the wave vector of the cooling laser has a component 
along the motional axis (cos O + 0) the integral of @ over the energies converges. Then, the value of F is found 
from J°° dEPo(E) = 1, yielding: 

F =m- (52) 

Using these results, the steady state energy has the form: 



(E) 



r°° n 
I ^W> = jci 



a + cos 2 0p / 7 2\S\\ 



4 cos 2 O \2\5\ 7 



Eq. ( p3| ) contains the dependence of the cooling limit on the angle 0o between the cooling laser beam and the 
direction of the motion. The final energy is minimal for cos0o = 1, i.e. when the laser propagates parallel to the trap 
axis, and it diverges for cos0o = 0, when the laser is orthogonal to the trap axis and there is thus no laser cooling. 

The minimum of the final energy vs. detuning is reached for S — —j/2, as in the case of one ion (see for example 0). 



Inserting into (53) N = 1, a = 1/3 (which corresponds to spatially isotropic spontaneous emission), and cos0 o = 1, 
we find the same result as Javanainen and Stenholm in their semiclassical expansion for one ion ||. Hence, the final 
energy for N ions is N times the steady state energy achieved by Doppler cooling of one ion. This general result 
has also been found earlier in special cases, such as a Coulomb cluster in the Lamb-Dicke regime p2| , and a two- ion 
crystal treated in Ref. |23] by extending the method of . 
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2. Time Evolution 



The steady state solution suggests to use the following ansatz for solving Eq. (|4^): 

P(E, t) = F{t)E N ^ exp ("^) . (54) 

where i/(£) is a positive function of time and F(t) is a normalization factor, such that at any instant t the relation 
J °° dEF(t)P(E, t) = 1 is fulfilled, i.e. 

F ® = rtW' (55) 

This ansatz corresponds to assuming that the distribution P(E,t) is always thermal with average energy 

(E(t)) = NU(t). (56) 
Substituting into Eq. (E9) and using ( [46| ) we get a differential equation for U : 

^i^-^-if C + iV#-*V (57, 



drC/VC 7 / ^ V U J \ U 
Since this equation must be fulfilled for all values of E, we find 

dU 1 



which has the solution: 

; .2 a_ r ' 



, . rr ^^cos^oLYO) \ 1 
U(t) = U cxp f jv * J + c ' (59) 

where iVC/ = (-E(O)) is the initial energy. For <->oowe recover from Eq. ( |59|) the steady state solution (|5l|). The 
rate at which the steady state is reached is: 



2lu r cos 2 9 \L'{0)\ _ M 8lo r cos 2 d Q n 2 j\5\ 
~N ~ ~N {AS 2 + 7 2 ) 2 



iv± uw H wo vote n u l /cnx 
1 c ° o1 — ^7 — ~ 7TB. i 2V2 ' 



where we have used (^J). Thus the cooling rate increases linearly with the number of driven ions M. It is largest if 
all ions are driven, as expected. 

Since the results derived so far agree precisely with those found earlier in specific cases, the general procedure which 
lead us from the quantum mechanical equations to the rate equation for the energy can be considered the common 
basis which underlies and unifies these earlier treatments. Furthermore, we have shown that the final energy of an 
iV-ion crystal is N times that of a single motional mode, while the cooling rate is M times that for the one-dimensional 
motion of a single ion, M being the number of ions driven by the laser. 

A final remark should be made on the assumptions leading to (|i^). In deriving the Fokker-Planck equation we 
have assumed ergodicity and Eq. p^), i.e. that p(E) varies negligibly on the recoil energy scale. The latter condition 
corresponds to the second-order expansion in fi of Ref. || . In that work the derivation of a Fokker-Planck equation 
was based on the limit of ovcrdamped oscillation, 7 ^> u, in order to adiabatically eliminate the excited state from 
the equations. Our derivation, however, does not necessarily imply this limit. Only in the case of N = 1 ion, which is 
not the main focus of our study, we must have AE 3> hv in order to fulfill the condition of a large density of states, 
D{E) 1, and thus for this special case the overdamped oscillator limit is a requirement for the validity of 



B. Lamb-Dicke regime 

When the atomic motion is well localized with respect to the laser wavelength (Lamb-Dicke regime), the Franck- 
Condon coefficients in (^|) can be approximated by their first order expansion in the parameter ((k ■ x) 2 ). For a single 
ion excited below saturation, the dynamics are described by a rate equation for the populations of the states with 
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vibrational number n, which can be analytically solved In this form, since n is proportional to the mechanical 

energy of the ion, the equation for cooling is an equation for the energy. For many ions, a set of rate equations can 
also be derived which describe cooling of each mode and which are decoupled, since a simultaneous change of the 
energy of two or more modes by scattering of a photon is of higher order in the Lamb-Dicke parameter (and thus 
takes place on a longer time scale) pO|Jl7| . 

Yet we show in this section that under some conditions we can derive a single equation for the total energy of many 
ions. Given the analytical simplicity of the model, this example is also instructive in order to see the limits of validity 
of the ergodic equation. 

Let us consider Eqs. in the Lamb-Dicke regime |24|. Here the reduction of this set of equations into the rate 
equations (^3|) is possible without assumption (^2|), since the terms (p^), are negligible either because they are 
rapidly oscillating or because their coupling to the populations is of higher order in the Lamb-Dicke parameter [ pl)[ . 
In first order in the parameter u^/vp, with [3 — 1, —,N, the modulus square of the Franck-Condon coefficient has the 
form: 



N 



N 



Kn|e^|k)| 2 » ]j5 na , ka (2n + l) 



(61) 



a=l 
N N 



3=1 



P=la=l,a^3 

Thus, given the initial state of motion |n) = |m, ...,njv), scattering of a photon involves three possible transitions: 
(i) the so-called carrier transition, of zero order in the Lamb-Dicke parameter, with final state |k) = |n); (ii) the 
red-sideband transitions, where |k) — \rii,...,np — l,...,njv), and (hi) the blue-sideband transitions, where |k) = 
| Til, np + 1, npf). The cases (ii) and (hi) are of second order in the parameter rjp. 
When we substitute (|6l|) into Eq. (23), we obtain: 



A' 



dt 



(n|p|n) = -(nlplnJicjRCos 2 ^^ [(rip + l)L(vp)/vp + npL(-vp)/vp] 



(62) 



p=i 



N 



(n\p\n)u; R J2 a O/3 + l)L(p)/v/3 



8=1 



N 



J2( n + 1 p\p\ 11 + V-^Ka + 1) [cos 2 e Q L{-vp) + aL(0)] 

3=1 9 



N 



+ l/8|p|n- lp)—np [cos 2 9 Q L{vp) + aL( 

3=1 9 



where have used Y] j rjj = ujR/fp, and the vector lp is defined in the ./V-dimensional Hilbert space through (lp) a = 
S a ,p for a = 1, N. 

Eq. (62) is linear in the vibrational numbers np of the single modes, which are weighted by the factors u)nL(yp)/vp. 
Only when these weights are of approximately equal magnitude, a hypothesis similar to (|27| ) , which allows summation 
over states of equal energy, can be applied. This is true when ^ « 7 for all modes (3 because in this limit the 
Lorentzian atomic resonance curve varies very slowly over all red (blue) sidebands of the modes. Then, by expanding 
the terms in ( |62] ) at the second order in the parameter ^3/7, we obtain a Fokker-Planck equation of the same form 
as (||). In that sense, the condition vp <C 7 can be considered a requirement for deriving an ergodic equation in the 
Lamb-Dicke regime. 

Alternatively, in [ fL^ a single equation for the total motional energy in the Lamb-Dicke regime has been justified 
by assuming that all modes thermalize on a faster time-scale and thus imposing 



IE 

{np) = NmTp 



(63) 



Application of (|6^) to ( p2| ) yields an equation for the total mechanical energy which can be solved analytically. In our 
treatment above, in the case where we find a Fokker-Planck equation, condition (|6^) is also fulfilled, but there it is a 
consequence of the average coupling among energy shells due to light scattering, rather than an extra assumption. 
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C. Discussion 



The Fokker- Planck equation (^9() has been derived in two different ways: in the semiclassical case by starting from 
the ergodic equation (|3l|), and in the Lamb-Dicke regime by starting from equation (p3|). Both derivations rely on 
the limit 7 3> v,ojr, while in the Lamb-Dicke regime there is the additional constraint that cur <C v. 

Let us now summarize the solutions of (^9|). The form of the solution is the same as the one obtained for cooling 
of one ion, whereby here the number N of ions appears in the steady-state energy (^) as a scaling factor, and the 
cooling rate (^) scales as the ratio M /N where M is the number of driven ions. This result supports the conclusion 
drawn in that the cluster is cooled like a single ion. It is instructive, at this regard, to compare two particular 
cases which exhibit analogies: The axial motion of two ions in a trap with collective modes at frequencies Ui, ^2, and 
the motion of a single ion in a two-dimensional harmonic oscillator potential with the same frequencies, v\ and z/2, 
on the two axes. Neglecting for the moment the different masses, the mechanical Hamiltonians of the two systems 
are equivalent. The laser-ion interaction terms are also equivalent, provided that one of the two ions in the chain is 
driven, and that the Lamb-Dicke parameter for each mode of the chain is the same as the Lamb-Dicke parameter for 
each axis of the trap of the single ion. This latter condition can be fulfilled by choosing, in the single-ion case, the 
proper angle of propagation of the laser in the two-dimensional plane. 

The main difference between the two cases arises in the spontaneous emission: in the two-dimensional one-ion case, 
the photon is scattered at a random angle in the plane and the ratio between the projections of the recoil energy 
on each axis can vary. In contrast, in the two-ion case, the scattering angle is the same for both modes, hence the 
average energy transfer to the two modes has a fixed ratio in any scattering event. This difference appears as a the 
geometrical factor in the dynamical equations, as well as in the expression for the final energy (|5^) [p^| . From that 
consideration we expect that the generalization of our treatment to a 3-dimensional Coulomb crystal will change the 
one-dimensional results only by numerical factors representing the different geometry of the problem. 



V. EVALUATION OF THE SEMICLASSICAL ERGODIC EQUATION 

So far we have assumed a grid of energy AE, which defines the coarse-graining on which the relevant system properties 
for cooling are defined. The grid AE has been chosen according to the condition h-f ^s> AE, and an equation for the 
energy has been derived for the situation that the number of states at energy E is D(E) 1. Thus, AE represents 
a limited resolution which allows us to average over many quantum states. We now discuss the limit AE 3> Tivjq- 
This case is expected to correspond to the classical limit, since the resolution is such that all details of the quantum 
spectrum are averaged out. We will derive the classical limit starting from the average Franck-Condon coefficients 
QW(E,E') in (§§), from which we obtain an explicit form for the ergodic equation (|3l|). 
Let us consider the sum in Eq. ( |26| ) . Using the property of the trace, we can write: 

£ £'£' ^ lkZ3 i n > ^ elkZj i k > = £ £'£' Tr I i k > M e ~ lkZ3 i n > ^ kZi ) 

j k n j k n 

= C(E k ,E n ), (64) 

thus defining C(E,E') = QW (E, E')D(E)D(E'). The summation over the states belonging to the shell of energy E 
can be re-expressed as: 



£'|k)(k|=/ d£5(£-JT mec )£|k)<k|, (65) 

k k 

and Eq. (|64|) aquires the form: 

C(E*, E n ) = 4 £ Tr If dE f dE ' d ( E - H mC c)e lkz >6(E' - H mcc )e~ tkz > \ . 



Using the Fourier transform of the 5-function, we get: 



dES(E - H mec ) = I dE^- r d ™*<*-*-W* = -1 r dT Sin(AErA) e^-^ 

lF{E k ) Jr(E k ) 27 ™ J -00 2n J-00 T 



r/h 



(66) 



(67) 
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such that Eq. 



can be rewritten as: 



(68) 



where we have defined: 



Ay (r, r') = T I { e i(E-6 mc )T/n e ikz ie i(E'-& n 



o)r'/n. 



1 



(69) 



By using coherent states ja/j) to calculate the trace, using (|4|), and the properties of coherent states [^lj, Eq. (|6 
takes the form: 



JY 



A i (r,r') = n 

0=1 



d 2 ag 



r 27V 



exp 



-if(r + r') 



exp 



-|a^| 2 e-^ r '(l-e-^ r ) 



• exp 



M 2 + + «7?Ke-^ - a/3 ) ) (1 - e"^' ) 



(70) 



In Eq. (pq), the integrand is the product of Aj(t,t') and the window function sin(x)/x of width h/AE. According 
to Eq. (J70|), Aj(t,t') depends on t,t' only through e lUf3T and e Il/ ' 3T , and thus r is scaled by the mode frequencies. 



A. Derivation of the classical limit 



We now consider the limit E, AE ^> Twp for all (3. This corresponds to the situation (np) ^> 1, e.g. before the cooling 
limit is reached, and in general to the regime 7 » i/^. In this limit we can expand the exponentials in (|7^) in the 
parameter hvj AE. Up to first order, Eq. (|6"8|) is equivalent to: 



C(E,E>) 



(¥) 

1 

2^ 



2 N 

n 

18=1 



d 2 ^ 



r 2JV 



- r dr 

2tt 



siii(Agr/fi) i(g/fi _ yw(|afl| 2 +1/2))T 



,sm(AET , /K) , (E 



dr' 

A£r'/ft 



(B7r l -^(|a^| 2 + l/2-27 ? ' 3 Im(a3)-r ; ' 32 ))T' 



(71) 



This result can be interpreted physically, considering that r is the Fourier conjugate of the energy E/h and has the 
dimension of a time. Thus, Eq. ( |69| ) can be seen as the overlap between the state \a) and the state \a'), which 
corresponds to a displacement of followed by free evolution for a time r, then by a displacement back, and 
finally by free evolution for a time r' . Integral ( |68| ) sums over all intervals of times t, t', for all trajectories in the 
neighbourhood of the energy shells E at r and E' at r', weighted by the time window of resolution h/AE. For 
t = t' = the overlap is maximum, since the state is unchanged. Although the overlap shows a certain periodicity, 
e.g. for N = 1 and r, r' multiples of 2w/v, A(t,t) = A(Q,0), for AE fiv these recurrence times fall out of the 
window function interval. Therefore the only appreciable contribution to the integral comes from r, r' ~ 0. In other 
words, the classical limit corresponds to large energy uncertainty, such that only the short time evolution of the wave 
packet contributes appreciably to the integral, since only in this limit the trajectories are phase coherent. 

Eq. (fn]) can now be rewritten as an integral in the classical phase space using the definition of coherent states 



where q/3,pp are real numbers. Integrating over t, t' yields: 

dqi...dq N dp 1 ...dp N 

h N 



C(E,E') = / cLEi / dE 2 



(72) 



(73) 



N -2 



2m 



5 E, 



(pi - hk ) 
2m 



2 H ^2 

-V 

2m 



n 



J2^-V(qt,...,q N ) 

(3=2 
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Integrating in phase space [see Appendix C], and using the relation Je(E') — C(E, E 1 ) / (D(E')AE) we obtain: 



fE(E') = 



r(N)/y/4huj R E 

V^r(iv-i) 



1 - 



E' -E-%uji 
^JATiujrE 



2V-3/2 



(74) 



The function /b(-B') is real, and thus well-defined, on the interval of energies [E+%u) R — s/^iJrE, E+ TujJr + -y/ ATwjrE] , 
and it is normalized with respect to E' . In Fig. 3 we plot Eq. (|74|) for N = 1, 10, 100. For N = 1 it has the well- 
known form of the classical momentum distribution of a harmonic oscillator at a given energy E: In fact in this case 
E 1 = E + hojR + hkp/m, i.e. the probability for the oscillator to have final energy E 1 after the scattering is equal to 
the probability that the oscillator, at energy E, has momentum p when the scattering event occurs. For higher TV the 
distribution peaks more and more narrowly around E + Tko R . 
Finally, note that, from the definition of fsiE') in (|74[) and from (38): 



Q (k) (E,E') 



.u N T(N) 



^T{N - \){Ahu R EE') 



N- 



-[-{E'-Ef 



fc2 2 

H UJ R - 



2Hlu r (E' + E)] 



N-3/2 



(75) 



This function is symmetric in E and E' as in the quantum mechanical case (see (|26|)). In the limit in which this 
result holds, we eventually find the explicit form of the ergodic rate equation by inserting Eq. ([75]) into Eq. (|3l|). An 
explicit form of the rate equation can also be evaluated in the limit of N ^ 1 ions p5| . 
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FIG. 3. Plot of g(E)Q (k \E,E') as in Eq. @ as a function of E for N = 1, 10, 100. E,E' are in arbitrary units. Here, 
£7=1, Tiujr — 0.5. The dashed line indicates the location of the center of the distribution, E 1 = E + Kuir. 



VI. CONCLUSIONS 



We have studied Doppler cooling of a Coulomb crystal: Starting from the full Master Equation, in the low saturation 
regime and by consecutive steps of averaging we have derived a rate equation for the total energy of the crystal, 
reducing dramatically the number of degrees of freedom and thus simplifying considerably the complexity of the 
problem. The equation is defined on a coarse-grained energy scale with grid AE such that AE <C Try, where 7 is 
the linewidth of the electronic transition resonant with laser light. Its derivation is based on the quasi-continuum 
characteristic of the spectrum of motional energies E on AE, and on the assumption that the coupling between states 
belonging to different energy shells is a smooth function of E. 

Starting from the general form of the equation, we have studied the scmiclassical limit and the Lamb-Dicke limit, 
and in both cases we have found a Fokker-Planck equation describing cooling of a Coulomb crystal of N ions. The 
general solution agrees with the results of previous treatments, which were developed in perturbation theory for these 
limiting cases ^,^,^,^2|. As observed in [Q, the dynamics of cooling of a Coulomb crystal can be scaled to the one 
of a single ion. 

In the semiclassical limit we have derived the explicit form of the rate equation for the mechanical energy, by 
calculating from the full quantum mechanical expression the classical probability of scattering between motional 
states at different energies. An explicit form of the equation can also be derived in the quantum limit for the case of 
N ^> 1 ions. This derivation will be presented in future work Ba). 
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Our results, with marginal changes, can be applied to a Coulomb crystal in three dimensions. In that case, the 
dimensionality enters into the number of modes, which for a crystal of N ions is 37V, and into the spatial distribution 
of the scattered photons, thus affecting the coefficients of the energy rate equation. From the comparison between 
the cooling problem for a one-dimensional crystal and for one ion in more dimensions, we expect the final energy in 
3D to differ from the one-dimensional result only by a geometrical factor. 

To conclude, we would like to remark on the generality of the treatment. Other incoherent processes in physical 
systems can be studied in an analogous way, provided that the rate determining the dynamics of interest can be singled 
out, and that on the corresponding energy scale the spectrum of energy levels is characterized by a quasi-continuum. 
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APPENDIX A: MOMENTS DERIVATION 



One ion. Let us consider a one-dimensional harmonic oscillator of frequency v and number state |n). Considering 
\(n\ e~>cp(ikx)\k}\ 2 as a distribution over the final states of energy E^, given the initial state \n) with energy E n , 
the first and the second moments of the distribution are: 

(£ k -£„> k =^^(fc-n)|(n|e^ at+a )|fc)| 2 , (Al) 

k 

((£ k - £ n ) 2 ) k = ftV ]T(fc - ,^ 2 |He"'( at + a )|£ ; }| 2 . (A2) 



Using the property a^a\k) = k\k) and the closure relation of the states we can contract the sum over k in (Al), 

(|A|), using: 

/HHe"' (at+Q) |fc)| 2 = (n|e~"'( at+a ) (a) a) m e"^ +a ^\n) . (A3) 

Expression (A3) can be further simplified by using the commutation properties of the bosonic operators, and we 
obtain: 

(£ k - £ n )k = Tiujr, (A4) 
((£ k -£ n ) 2 ) k =2htu R £ n + h 2 Lj 2 R , (A5) 

where we have used rj 2 = wr/v. 

N ions. Let us now take an A-ion chain. Using the definitions of Section II, the first moment of the distribution is: 

m-Snh^ =^(f k -f n )|(k|e^|n)| 2 

k 

2 



(A6) 



where the subscript j refers to the driven ion. Analogously, the second moment has the form: 



(<(£k - £n) 2 )k) . = (^-) v^v a {2n a + 1) (b«f + h 2 u; 2 R . (A7) 



Averaging (A7) over the ions of the chain we find: 

' 3 
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APPENDIX B: EVALUATION OF THE SEMICLASSICAL TRACE 



Expression ( |73| ) can be rewritten as: 
F = / &Ei / d£ 2 / dptS (e 2 -Ek- hcj R ~ —ft") / dgr-dgjv ^'^ ^ ( ^ - V %L - V{q u ..,q N ) 



8=1 

(131) 



We define: 



I(E, Pl ) = / dgx..^ dp2 "^ Py ^))' ( B2 ) 

and move to the coordinates qg where V is a diagonal quadratic form. By introducing the set of rescaled variables 



Qg — ^mi'p/2qa, Pg — y/l/2mpg, integral ( jBg ) is the measure of the surface of a unitary hypersphere in 2N — 1 
dimensions. Integrating, we obtain: 

_ 2 N {E-pl/2m) N -^ 

mpi) = ^ V1 .., N r(iv-i/2)" (B3) 



Substituting now this expression into (Bl) and integrating over p\ we get: 

e=[ m , 1 - J d S2 -_l ( Ei _iE 2 ~n. R -E l? ^ 

J H E) T(N - 1/2) } nBI) 2 hk^A N Vl ...v N \ ihcuR 

dE ^TN 1 ^ / dE 2 f El (E 2 ), (B4) 



where 



F(E) h Vi...V N T(N) J3={E') 



0Fr(7V - 1/2) y/Ahio R E V 4?iu; fl £; 



It can be easily verified that (B5) satisfies the relations fl38|) and fl39|-|41 
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